Cellular Models for River Networks 



O 

o 
o 

(N 
> 

o 



Guido CaldarcUi 

INFM, Sezione di Romal, Dip. Fisica, Universita di Roma "La Sapienza" , P.le A. Moro 2, 00185 Roma Italy 

A cellular model introduced for the evolution of the fluvial landscape is revisited using extensive 
numerical and scaling analyses. The basic network shapes and their recurrence especially in the 
aggregation structure are then addressed. The roles of boundary and initial conditions are carefully 
analyzed as well as the key effect of quenched disorder embedded in random pinning of the landscape 
surface. It is found that the above features strongly affect the scaling behavior of key morphological 
quantities. In particular, we conclude that randomly pinned regions (whose structural disorder bears 
much physical meaning mimicking uneven landscape-forming rainfall events, geological diversity or 
heterogeneity in surficial properties like vegetation, soil cover or type) play a key role for the robust 
emergence of aggregation patterns bearing much resemblance to real river networks. 
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I. INTRODUCTION 

Through experimental studies, it has become evident 
in the past few years that the geometrical and topolog- 
ical structures of river basins are characterized by the 
absence of a single well-defined length scale. This is re- 
flected in the appearance of power laws in the distribu- 
tion of several quantities, chiefly total contributing area 
at a point j^] and stream lengths |^-^, and by the clear 
experimental assessment of scaling properties (yielding 
either self-similarity or self-affinity pi) for many geomet- 
rical descriptors of the river basin . The discovery 
of the general underlying mechanisms yielding scale-free 
features is the present theoretical challenge. 

The network associated to a given natural terrain per- 
taining to a river basin can be experimentally analyzed 
by using the so-called Digital Elevation Map (DEM) tech- 
nique P,p|-p^ which allows to determine the average 
height of areas (pixels) of the order 10~^Km^. Thus a flu- 
vial basin is represented in a objective manner often over 
4 log scales of linear size. Lower bounds are imposed by 
channel initiation processes at 0(10—100) m. Crossovers 
of geological nature provide altered aggregation processes 
and thus an upper cutoff, usually beyond scales of 0(10^ 
- 10^) m. Thus the observational evidence yields a much 
reliable framework over many scales for comparison with 
dynamical models aimed at the origin of scale-free fea- 
tures. 

Much interest has been recently attracted by land- 
scape evolution models. Chief among those are the de- 
tailed deterministic models which address the description 
of the detailed dynamics acting on the landscapes 
The reductionist approach, where a precise description of 
the details of the dynamics is sought, is successful, and 
much interesting, in the pursuit of the description and 
the classiflcation of landforms. Nevertheless, as standard 
in critical phenomena, the mechanism producing scale- 
free structures is expected to depend only on a few key 
features common to all the networks rather than on the 
details of the particular system under study. Hence in 
this work, centered on the dynamic origin of fractal river 
networks, we follow a nonreductionist approach based on 



the simplest possible, parameter-free models capable of 
allowing the emergence of complexity. 

A flowrate unit is associated with each pixel and the 
flow contributing to any pixel follows the steepest descent 
path through drainage directions whose collection deflnes 
the planar structure under consideration. The resulting 
network is therefore the two-dimensional projection of 
the three-dimensional tree-like structure of the steepest 
descent paths draining a given basin. The planar pat- 
terns of network aggregation are obtained by employing 
the cellular model for the evolution of a fluvial landscape 
originally introduced by ||l^ and further studied in jl^] 
It is aimed at describing in a crucially simple manner 
the sole fluvial component of landscape evolution. Al- 
though such component must be coupled to other - chiefly 
hillslope - transport processes to yield a comprehensive 
dynamical description pd|JI^ , p7| -pJt , it rules the planar 
imprinting of the network. Hence the detailed study of 
the model is deemed signiflcant. 

Starting from a three dimensional landscape, evolu- 
tion occurs according to a threshold dynamics similar to 
the one proposed in self organized critical (SOC) mod- 
els [ pO| . The main idea of the erosion dynamics is that 
whenever the local shear stress exceeds a given threshold, 
erosion starts an 'avalanche' and a related rearrangement 
of the network patterns takes place. The model of self- 
organizing fluvial structures may be seen as a modifica- 
tion of the sandpile model developed as a paradigm of 
the dynamics of open, dissipative systems with many de- 
grees of freedom. It may be thought of as belonging to the 
set of models in which the threshold for activity, say Tc, 
rather than being a constant value, depends on non-local 
properties of the self-organizing structure. In the fluvial 
case, the non-local character of the threshold value fol- 
lows from the fact that the threshold at the arbitrary site 
equals a shear stress, i.e. r cx Vhy/a, where h is the local 
landscape elevation and a is total contributing area sur- 
rogating total flow collected from a distributed rainfall 
event. As such, the exceedence of Tc depends not only 
on local conditions (i.e., a critical value of V/i), but also 
on non-local conditions defined by the contributing area 
a computed through drainage directions, i.e. it depends 
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on the entire state of the system which is self-organizing. 
Notice that the physical rationale for the nonlocal depen- 
dence lies in the fact that the system is open, i.e. injected 
from outside, allowing flow rates to be proportional to to- 
tal contributing drainage area. The long-range nature of 
the threshold dynamics tends to hide from the observer 
the temporal fluctuations which take place in the evolu- 
tionary time scale. In this sense the above model was 
classified Q as one of spatial self-organized criticality. 

Whether or not river self-organization qualifies as a 
more general framework of self-organized criticality re- 
mains to be seen. If SOC must necessarily refer to the 
occurrence of a critical state in the sense of critical phe- 
nomena, where a small local perturbation can cause a 
significant change in the configuration of the whole sys- 
tem and thus the system shows both spatial and tempo- 
ral scaling, then the time dynamics should be specifically 
considered. One way to do this is through the oscillation 
of the threshold in time, i.e. Tc{t), simulating climatic 
fluctuations (see |^), through which indeed temporal 
evolution appear, or through perturbations of random 
location and strength in the evolution of the landscape. 
This is also true, as we will discuss later, if the landscape- 
forming rainfall events are described as nonuniform in 
space, leading to patches of activity randomly scattered 
spatially (in such a case the outflow response of the sys- 
tem becomes a 1// signal). However, regardless of any 
additional features, we believe that the central scope of 
SOC is the dynamic explanation of the growth of frac- 
tal structures of the type appearing in nature, i.e. the 
physics of fractals. As such we feel that our classification 
of the model as a particular case of SOC is a suitable 
one regardless of the description of the embedded tem- 
poral activity because the system always reaches a frac- 
tal state. Moreover, questioning on this basis |p3|] t he 
self-organized critical nature of the model by Ref. |15[| is 
irrelevant because it has been shown on thermodynam- 
ics grounds that scaling properties of energy and entropy 
yield limit states which, depending on the constraints, 
are temporally frozen or active 1^. Furthermore, ref. 
shows that optimal states like the ones dynamically ac- 
cessed by the above model may exist in temporally active 
states precisely at the edge of a chaotic behaviour. 

A question, indeed more interesting than the semantics 
of SOC, is whether the constraints in the model may be 
relaxed to produce a 'hot' fluvial landscape more closely 
resembling an ordinary sandpile. This question is ad- 
dressed in H] and in more dctatil in 

Our main goal is twofold. On one hand we will extend 
previous investigations both in accuracy and in statistics 
by performing simulations at much larger scales. On the 
other hand, we will consider important issues such as the 
effect of the boundaries and of the initial conditions bear- 
ing much significance on geologic influences. In particular 
we will show that both the aforementioned effects play 
an important role in the results previously obtained. We 
also study the effects of disorder, say through the pres- 
ence of small, uncorrelated inhomogeneities in the initial 



conditions, in particular with regard to the robustness to 
single/multiple outlet arrangements. 

The paper is organized as follows. In the first section 1 
the model is recalled. Section 2 presents the results with 
emphasis on scaling analyses, while the following section 
focuses on the important effects of quenched pinning on 
the system. A set of conclusions closes then the paper. 

II. THE MODEL 

We consider a lattice model of a real landscape. Let 
be the height of the landscape associated to every site 
X of a square lattice of size L x L. The lattice is tilted 
at an angle 9 with respect to a given axis to mimic the 
effects of gravity. Two possibilities will be analyzed: 

1 All the sites on the lowest side (kept at height 
h — 0) are possible outlets (i.e. the multiple out- 
let arrangement) of an ensemble of rivers which are 
competing to drain the whole L x L basin; 

2 Only one site is kept at /i = and it is the outlet 
of a single river in the L x L basin. 

In addition, for both the above cases, two types of 
initial conditions will be considered: (a) a regular initial 
landscape, e.g. flat, and (b) an irregular surface obtained 
by superposing to a smooth sloping surface a suitable 
noise. 

Each site collects an unit amount of water from a dis- 
tributed injection (here a constant rainfall rate as in the 
original approach) in addition to the flow which drains 
into it from the upstream sites. A unit of water mass is 
assigned to each pixel of drainage area so that the total 
area drained into a site is also a measure of the total wa- 
ter mass collected at that site. From each site water flows 
to one of the eight sites, four nearest neighbors and four 
next nearest neighbors, having the lowest height (i.e. the 
steepest descent path). We shall indicate all these eight 
neighbouring sites as nearest-neighbours (nn). This con- 
struction allows the assignment of drainage directions to 
an arbitrary landscape. The drained area ax is associated 
to each site x according to the equation 

where the sum runs over the subset y{x) G nn{x) of 
neighbor sites whose area is actually drained by x. The 
second term in eq.(|l|) represents the uniform injection. 

Also, the (up)stream-length from site x to the source 
is computed according to the following procedure. At 
a given site x the areas of all nn{x) of that sites are 
checked. Following the ordinary meaning of downstream 
and upstream sites (i.e. downstream is the site one finds 
following the river to the outlet, upstream is the site fol- 
lowing which one reach the source from which the largest 
incoming river enters the site.) Following Ref. the nn 
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with largest value leads to the outlet and is defined to 
be a downstream site. The nn with the second largest 
value indicates the longest path toward the source and 
is defined to be the upstream site. The sum of all the 
upstream sites from site x to the source is Ix ■ The down- 
stream length could be defined through an analogou pro- 
cedure. Experimental measures are available for both ax 
and §. 

The time evolution of the model follows the following 
steps: 

1 The shear stress Tx acting at every site is computed 
according to [|l5| 

Tx = Ahxy/Chi (2) 

where Ahx is the local gradient along the drainage 
direction. 

2 If the shear stress at a site exceeds a threshold 
value, Tc, then the corresponding height hx is re- 
duced (i.e. by erosion) in order to decrease the 
local gradient. The shear stress is set just at the 
threshold value. This produces a rearrangement of 
the network followed by a reupdating of the whole 
pattern as in step 1. 

3 When all sites have shear stress below threshold 
the system is in a dynamically steady state. Since 
this situation is not necessarily the most stable, a 
perturbation is applied to the network with the aim 
of increasing the stability of a new steady state. A 
site is thus chosen at random and its height is in- 
creased in such a way that no lakes, i.e. sites whose 
height is lower than that of their eight neighbours, 
are formed. Steps 1 and 2 then follow as before. 

After a suitable number of the perturbations (step 3), 
the system reaches a steady state which is unsensitive 
to further perturbations and where all statistics of the 
networks are stable. This resulting state is scale-free, 
i.e. it is characterized by power-law distributions of the 
physical quantities of interest. 



Model B An inclined plane with all sites at the bot- 
tom of the plane allowed to be possible outlets. 
This choice was selected with the aim of investi- 
gating the differences arising when arranging the 
boundary conditions with multiple outlets versus 
single outlet. The former allows for competition 
for drainage area among rivers. 

Model C,D The two previously considered situations 
with the addition of a random, uncorrelated noise 
(whose strenght - i.e. variance - is less than 10% of 
the average height). That is, on top of the height 
computed according to the rules of model A and 
B respectively (a comb-like lattice and an inclined 
plane), we added a random Dh that is extracted in 
the interval [- < h > /lO, < h> /lO] where < h> 
represents the mean altitude of the landscape. 

An average over a few (up to five) configurations was 
taken. This choice, especially when coupled to large sizes 
of the system, proves sufficient for a statistical descrip- 
tions sought in view of the self-averaging nature of the 
random perturbation. 

In Figure 1 typical landscapes sculpted by the above 
dynamical process and the corresponding networks 
drawn through the steepest descent construction are 
shown for models A,B. The same picture for models C,D 
is shown in Figure 2. 

Two features can be grasped from these pictures. First, 
in the results of both model A and B there is a strong 
memory of the initial configuration despite the fact that 
the dynamics of the erosion process was somewhat ex- 
pected to be sufficiently strong to soon lose the imprint- 
ing of its initial condition. Secondly, the single outlet 
restriction imposed in model A appears to be a severe 
constraint because it increasingly affects the wandering 
of the main river towards the lowest part of the basin. 
Our results suggest that this is indeed the case for flat 
initial conditions (A and B) while for noisy initial condi- 
tions boundary effects are of lesser importance. 

B. Area and length exponents 



III. RESULTS 

A. Landscape evolutions 

Our numerical calculations were carried out on a 
bidimensional square lattice (where each site has eight 
nearest-neighbours) for sizes up to L = 200 with reflect- 
ing boundary conditions in the direction transversal to 
the flow and open boundary condition in the parallel one. 
We considered the following initial conditions: 

Model A A comb-like structure with a single outlet 
This was the situation originally studied in [ p^ and 
our results are in agreement with theirs. 



Let us define P{a, L) and n(Z, L) as the exceedence (cu- 
mulative) probability distributions of the drainage area 
a and stream length I respectively arising in a domain of 
linear size L. The following scaling forms are expected 
to hold §: 

P[a,L)^a'-^F{j^) (3) 

n(/,i:) = ?i-^G(-^) (4) 

Here H is the Hurst exponent and di is the stream- 
length (or chemical distance) fractal exponent. 
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As it was already noted Q , for self-afRne river networks 
{H < 1, d/ = 1) the scaling relations relate all exponents 
in terms of H. For self-similar river networks {H = 1, 
di > 1), the same happens in terms of di. 

Another important indicator of basin morphology is 
the relation between the mean total contributing area a 
and the length of the main stream /max oc L'^' [p|,p4t, 
which is commonly known as Hack's law lE2|: 



ji/h 



(5) 



The related exponent has been studied in all simulations. 
A summary of the scaling relations between the various 
exponents involved is reported in Table |[ 

Experimental values of r and 7 are available from ear- 
lier analyses of DTMs from basins of different size, ge- 
ology, exposed lithology, climate and vegetation |p],p^. 
It was observed that, while a majority of basins tend 
to seemingly universal values r = 1.43 ± 0.02 and 7 — 
1.8±0.1, exceptions are observed where altered values are 
observed although always in a concerted manner. Since it 
was suggested Q that scaling laws for river networks are 
related, e.g. 7 = 1-|-(t— it was concluded there that 
no universal exponents are expected in nature. Rather, 
the roles of geology and tectonics concert a coordinated 
scaling structure which strives for fractality yet adapted 
to its geological environment. The results of the model 
described here, revisited in the above light, conform to 
this view. 

The results for the four models A,B,C and D for the 
area distributions are shown in Figure 3. It is appar- 
ent that, due to the pathological initial conditions, the 
scaling behaviour for models A and B is somewhat more 
noisy than for models C and D. Figure 4 contains the 
collapse plot for all the cases. Figure 5 shows the stream- 
length distribution for the four models. For this picture 
the same remarks of Figure 3 apply. In Figure 6 we show 
the collapse plot corresponding for the stream-length dis- 
tributions. 

A summary of the scaling exponents obtained is in- 
cluded in Table |l^ where we observe a consistent picture 
of related scaling exponents as theoretically expected - 
see Table |. 



C. Energy Dissipation and Optimal Channel 
Networks 

During the evolution of the landscape we also moni- 
tored the change in total energy dissipation of the system, 
defined as E = a°-^ (where x spans all sites of the 
lattice) ||25|-p7[| . The reason of this name come from the 
computation at any site of the gravitational energy lost 
by the falling of the water. In any point x one can expect 
a gravitational energy loss of the order of ax^hx where 
lS.hx represents the local gradient along the drainage di- 
rection. By using the observed scaling A/ij, cx ax'^'^ one 
obtains the above formula. The interest in this quantity 



comes from the fact that an extensively class of mod- 
els known as Optimal Channel Network (OCN) models, 
assume this quantity is minimized by natural landscape 
evolution. By using this principle OCN describe evo- 
lution from random spanning graphs to network more 
similar to the real ones. It is interesting to note that in 
this model where no hypothesis is made on E, we still 
observe an almost monotonical decrease of E associated 
with the dynamical evolutions, and a stabilization on dif- 
ferent plateaus of values of E. The actual figures for a 
sample 30 are as follows: E starts from an initial value 
of 7600 and decreases towards a plateau of 6800 — 6700 
where this monotonic decrease becomes slower (1% de- 
crease in 50000 steps). This behavior, also observed in 
other models [^8| , bears important consequences in the 
light of the suggested connection of fractality and opti- 
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IV. GEOLOGICAL CONSTRAINTS AND 
QUENCHED RANDOM PINNING 

This section presents a detailed study on the effects on 
landcape evolution of quenched randomness, simulated 
by a random choice of sites unable to evolve regardless 
of the threshold value developed therein. It is found that 
this form of disorder tends to favour aggregation pat- 
terns characterized by values of r = 1.43 ± 0.02 for both 
models, say, A and B (i.e. with single outlets or open 
boundary conditions). This results suggests that the ori- 
gin of the recurrent values observed in nature could be 
related to the ubiquity of heterogeneity in surface prop- 
erties characterizing locally the critical shear threshold. 

Within the river basin morphological and geological 
constraints play a definite role in the dynamical evolu- 
tion of landforms. The effects of quenched constraints, 
simulating any heterogeneity is the distribution of surfi- 
cial properties affecting erosion properties, is to favour 
some sites for the flowpaths, thus excluding other sites 
from the capture of the developing network which ulti- 
mately shapes the evolution process. In order to mimic 
such effects we analyzed the effect of a random pinning 
of a small region of the total surface (typically 5 — 9%) 
where the evolution is frozen, that is, the height is pinned 
to its initial value. We find that this constraint tends to 
favour aggregation even in the presence of random initial 
conditions. 

As regards the effect of the pinning, figure 7 (on the 
left) shows a sample whose dimension is 100 with mul- 
tiple outlet and random initial noise (model D). Figure 
7 (on the right) shows the same configuration (evolved 
from the same initial conditions) but with a 5% dilution 
pinning. It is evident that some of the smaller streams on 
the left have increased their size thus leading to a bigger 
aggregation. 

In this case we found t = 1.43±0.02. Moreover, all the 
other exponents verify the correct scalings predicted in 
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Table I. Purely for comparison purposes with the prevous 
results, we also report the plot of the area distribution 
in Figure 8 on the left. We also found that model D 
reproduces the same results, i.e. t = 1.43 ± 0.04 and 
7 = 1.60 ± 0.04. We are confident that this result at 
least for model C is quite robust with respect to changes 
of the pinning dilution. In the case of a 9% dilution for 
a smaller number of simulation we found quite similar 
results T = 1.44 ± 0.05 and a cumulative plot of P{a' > 
a, L) is shown on the right part of Fig. 8 p9|| . 

This result suggests that the origin of recurrent values 
observed in nature could indeed be related to the ubiquity 
of geological and morphological constraints in the surface 
properties locally characterizing the critical shear stress. 

V. CONCLUSIONS 

In this paper we revisited the model originally intro- 
duced in jisl l which we extended both in accuracy and 
goals. Specifically, we analyzed the stability of the uni- 
versality class of the original model with respect to the 
initial conditions and to the change from single to mul- 
tiple outlets. We found that if one starts with struc- 
tured initial conditions, critical exponents are sensible 
to a change from single to multiple outlets. On the 
other hand, upon starting from disordered initial con- 
ditions, we found critical exponents belonging to a new 
class which appears to be robust to the change from sin- 
gle to multiple outlet. Thus this simple model, under 
controlled conditions, yields somewhat different yet in- 
ternally consistent scale-free fluvial landforms depending 
on the dominant conditions affecting evolution. 

The above results conform to the experimental obser- 
vation p suggesting that the relevant scaling exponents 
for river networks are not universal. Rather, the frac- 
tal nature of river networks adjusts to the constraints 
imposed by the geological environment in a coordinated 
manner. It is interesting to observe that the final state 
of all simulation yields indeed fractal structures, as ob- 
served in nature, though characterized by different aggre- 
gation properties. The exponents characaterizing these 
different aggregates, nevertheless follow in a rather good 
agreement the scaling relations shown in Table I. 

We suggest that the lack of robustness in the value of 
the scaling exponents with respect to boundary and ini- 
tial conditions is related to the non-local character of the 
shear-based threshold, differently from what is observed 
in classical sandpile models of self-organized criticality 
|pO| . The intrinsic interest of the different aggregation 
properties of the steady states of the dynamics is related 
to their optimality with respect of total energy dissi- 
pation [ p^ . In fact, depending on external conditions, 
the same dynamical process may indeed get trapped in 
steady state configurations yielding local minima of the 
total energy dissipation functional, in what we may de- 
fine as a feasible optimality process ||^ . 



Finally, we have found that quenched disorder, mod- 
eled by random pinning, has a profound effect on the 
robustness of the resulting planar patterns by favour- 
ing aggregation and by locking the planar landforms into 
modes quite similar to the ones observed in nature. 

The remarkable success obtained by such a simple 
model in enlightnening some crucial features of the real 
basins is promising for a future success in a general char- 
acterization of the dynamics of fractal growth. 
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Exponent 


Self-Similar 


Self-Affine 


T 


2-di 


{1 + 2H)/{1 + H) 


7 


2/di 


l + H 


h 


di/2 


1/(1 + ^) 



TABLE I. Scaling relations: all the exponents can be 
determined in terms of di in the fractal case and H in the 
self-afBne case. 





Model A 


Model B 


Model C 


Model D 


T 


1.43 ±0.03 


1.50 ±0.03 


1.38 ±0.02 


1.38 ±0.02 


1 


1.60 ±0.05 


1.70 ± 0.05 


1.60 ± 0.02 


1.60 ± 0.02 


hth 


0.72 ±0.05 


0.71 ±0.05 


0.63 ± 0.05 


0.63 ± 0.05 


hex 


0.72 ±0.05 


0.69 ±0.02 


0.65 ±0.02 


0.65 ±0.02 



TABLE II. Data computed for the various computer sim- 
ulations. Since the scaling relationship h = one can also 
compute the theoretical value hex with the measured one h^x- 
The agreement of the two is rather good. 






FIG. 1. The final landscape (on the left) and the final 
network structure (on the right)for Model A and model B. 
Models A and B, start from a deterministic initial condition. 
Model A has single outlet, model B has multiple outlets. 




^1/ 




FIG. 2. The final landscape (on the left) and the final net- 
work structure (on the right)for Model C and model D. Mod- 
els C and D, start from a random initial condition. Model C 
has single outlet, model D has multiple outlets. 
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FIG. 3. Log-Log plot of the area cumulated distribution 
P{a,L) versus a for model A,B,C and D. The full line has 
a slope corresponding to r = 1.43, r — 1.50, r — 1.38 and 
r = 1.38 respectively. 
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FIG. 6. Scaling function Z^"^n(;,L) versus //L"' for mod- 
els A,B,C,D. The values used to obtain the collapse were to 
7 — 1.6, 7 = 1.7, 7 — 1.6 and 7 — 1.6 respectively. 
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FIG. 4. Scaling function a^~'^P{a,L) versus a/L^^^ for 
model C. The used values to obtain the collapse were 
T = 1.43, r = 1.5t = 1.38, and r = 1.38 for Model A,B,C,D 
respectively, = 0.6 in all the cases 




FIG. 5. Log-Log plot of the lenght cumulated distribution 
P{l,L) versus / for model A,B,C and D. The full line has a 
slope corresponding to 7 = 1.6, 7 = 1.7, 7 = 1.6 and 7 = 1.6 
respectively. 




FIG. 7. Comparison between the evolution of two identical 
initial configurations of model C with size L — 100 without 
(left) and with pinning (right). The pinning diluition was 5%. 




FIG. 8. Log-Log plot of the area distribution P{a, L) vs. a 
for model C with a 5% dilution (left) and with a 9% dilution 
(right). The slope of the lines corresponds respectively to 
r = 1.43 and r = 1.44. 
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